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Abstract 


In the center and south of the Mexican Republic, each year the hurricanes 
of the Caribbean Sea and the Pacific Ocean cause floods that lead to a 
wet season and that generally increase in magnitude and danger as the 
cyclone season progresses. Both conditions allow bivariate frequency 


analysis of their dates of occurrence and their maximum flows (Qm). In 
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this study, the bivariate distribution was formed based on the Gumbel- 
Hougaard Copula function, which satisfies the observed dependency 
condition (A$“) and which combines the von Mises distributions as 
marginal distributions for the dates of occurrence in the year and for the 
Qm a suitable probabilistic function. The exposed theory is applied to the 
annual floods recorded at the Guamúchil gauging station of Hydrological 
Region No. 10 (Sinaloa), Mexico, in the period from 1940 to 1971. The 
von Mises distribution is fitted via numerical optimization with the de 
Rosenbrock algorithm and the ideal distribution of the Qm turned out to 
be the Kappa. The graph of joint return periods of the AND type of 50, 
100 and 500 years was formed. In addition, conditional joint return 
periods of occurrence dates were estimated given that the Qm has the 
cited return periods. This allows estimates of the probability of 
exceedance of Qm in defined periods. The conclusions highlight the 
simplicity of these bivariate frequency analyses, by means of the Copula 
functions, and the practical importance of their predictions, according to 


the dates of occurrence. 


Keywords: Dates of occurrence, von Mises distribution, Copula functions, 
Kendall's tau ratio, joint empirical probabilities, dependency on the 


extreme right, joint and conditional return periods. 


Resumen 


En el centro y sur de la república mexicana cada año los huracanes del 
mar Caribe y del océano Pacífico originan crecientes que definen una 
estación húmeda, y que en general aumentan en magnitud y peligrosidad 


conforme transcurre la temporada de ciclones. Ambas condiciones 
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permiten el análisis de frecuencias bivariado de sus fechas de ocurrencia 
y sus gastos máximos (Qm). En este estudio, la distribución conjunta se 
formó con base en la función Cópula de Gumbel-Hougaard, que satisface 
la condición de dependencia (Af"“) observada y que combina como 
distribuciones marginales la de von Mises para las fechas de ocurrencia 
en el año y para los Qm una función probabilística idónea. La teoría 
expuesta se aplica a las crecientes anuales registradas en la estación de 
aforos Guamúchil de la Región Hidrológica No. 10 (Sinaloa), México, en el 
periodo de 1940 a 1971. La distribución de von Mises se ajusta vía 
optimización numérica con el algoritmo de Rosenbrock y la distribución 
idónea de los Qm fue la Kappa. Se formó la gráfica de periodos de retorno 
conjuntos de tipo AND de 50, 100 y 500 años. Además, se estimaron 
periodos de retorno conjuntos condicionales de fechas de ocurrencia, 
dado que el Om tiene los periodos de retorno citados. Lo anterior permite 
estimaciones de la probabilidad de excedencia del Qm en lapsos definidos. 
Las conclusiones destacan la simplicidad de estos análisis de frecuencias 
bivariados por medio de las funciones Cópula y la importancia práctica de 


sus predicciones, según las fechas de ocurrencia. 


Palabras clave: fechas de ocurrencia, distribución de von Mises, 
funciones Cópula, cociente tau de Kendall, probabilidades empíricas 
conjuntas, dependencia en el extremo derecho, periodos de retorno 


conjuntos y condicionales. 
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Introduction 


Generalities 


In general terms, the risk of flooding and its damages are a direct function 
of the water volume that flows through the channel which exceeds its limit 
or capacity, overflowing and covering the floodplains. However, the 
occurrence date of the event is as important as its magnitude and in some 
cases greater, when they occur out of season or in the wet season, 
because it takes the population off guard, causing greater damage 
(Khedun, Singh, 8 Byrd, 2019). 


Due to the above, the estimation of the probability of flood 
occurrence throughout the year, is vital for the elaboration of plans that 
have no hydraulic works of damage mitigation, which include the 
preparation for the event, with the purpose of reducing the exposure and 
vulnerability of the population, as well as optimizing the economic 
resources available for the emergency and accelerating post-event 
recovery (Durrans, Eiffe, Thomas, € Goranflo, 2003; Khedun et al., 
2019). 


In general, understanding the seasonal behavior of floods is vital in 
the planning and management of the river's hydraulic resources, both for 
agricultural and hydroelectric uses, as well as for navigation, recreational 
uses, and other activities associated with waterbodies. Therefore, 


knowing the relationship between the maximum flow and its date of 
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occurrence becomes extremely important, to the point of requiring its 


joint bivariate study (Chen, Guo, Yan, Liu, 8 Fang, 2010). 


Bivariate flood frequency analysis began formally at the beginning 
of this century (Yue 8 Rasmussen, 2002) and was generally based on 
maximum flow and of annual flood volume, using the bivariate Normal 
distribution and the Logistic model which accepts as probability 
distribution functions (FDP, by its acronym in Spanish) equal marginals to 
the extreme value distributions, the most common being the Gumbel and 
the GVE (Escalante-Sandoval 8, Reyes-Chávez, 2002; Aldama, Ramírez, 
Aparicio, Mejía-Zermeño, €: Ortega-Gil, 2006). 


Bivariate flood analysis is the simplest multivariate approach, and 
yet it involves five mathematical complications: (1) a bivariate FDP must 
be used; (2) its validation requires the estimation of bivariate empirical 
probabilities; (3) now there are joint and conditional probabilities; (4) a 
joint return period must be defined, for which there are infinite pairs of 
values of X and Y that can satisfy it and (5) the critical or design events 
must be selected among the aforementioned pairs (Ramírez-Orozco 8 
Aldama, 2000; Escalante-Sandoval 8, Reyes-Chávez, 2002; Volpi 8, Fiori, 
2012; Requena, Mediero, € Garrote, 2013). 


Currently, with the use of the mathematical tool known as 
"Copulas", the bivariate FDPs can be constructed with marginals of 
different types, given the fact that copula functions allow multivariate 
distributions to be represented from univariate or marginal FDPs, 
regardless of their shape or form type (Salvadori, De Michele, Kottegoda, 
8: Rosso, 2007; Meylan, Favre, 8 Musy, 2012; Genest €: Chebana, 2017; 
Zhang € Singh, 2019; Chowdhary 8 Singh, 2019). 
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By means of the Copula functions, the bivariate frequency analysis 
of the annual flood occurrence dates and their maximum flows is 
addressed. The former are represented by the von Mises distribution and 
the latter by a suitable PDF. 


Objectives 


The objectives of this study can be encompassed by the following six: (1) 
the directional statistics that represent the flood occurrence dates during 
the year are exposed; (2) the von Mises distribution (dvM, by its acronym 
in Spanish) and its adjustment via numerical optimization are described, 
which allows the probabilistic characterization of the dates of occurrence; 
(3) the basic characteristics of the Frank and Gumbel-Hougaard copula 
functions (FC, by its acronym in Spanish) are presented, including: 
Kendall's tau ratio, observed and FC dependence, estimation of joint 
empirical probabilities, and selection and ratification of the FC; (4) the 
selection and adoption of the ideal marginal FDP of the maximum annual 
flows is exposed; (5) the joint return periods of the OR, AND and 
conditional type are described and (6) as a numerical example, the 
application of the theory and procedures exposed to the 32 annual floods 
registered in the Guamuchil hydrometric station of the Hydrological 


Region No. 10 (Sinaloa), Mexico, is detailed. 
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Operative theory 


Circular data and directional indicators 


The dates of occurrence of the annual floods in Mexico, generally define 
a wet season that covers from June to October, which is the season of 
hurricane incidence that are generated in the Caribbean Sea and in the 
Pacific Ocean. Therefore, the dates of of the annual flood occurrence can 


be represented within a circle, which covers the 365 days of the year. 


There are several conventions or ways of working in the circle, to 
graph the data (Ramiírez-Orozco, Gutiérrez-López, € Ruiz-Silva, 2009). 
From now on, the Burn convention (Burn, 1997) will be used, due to its 
similarity with Cartesian quadrants. In such a scheme, the advance is 
counterclockwise, starting on the abscissa axis; for this reason, January 


1 and December 31 coincide in such beginning. 


Having several circular data drawn, it is possible to obtain its 
directional indicators, the most important are two, its mean direction (x) 
and its seasonality index (r). The first defines the central tendency of the 


data and, therefore, is the average occurrence date of the annual floods 


and the second quantifies the dispersion of such values (Campos-Aranda, 
2017). 
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Average address and seasonality index 


To estimate these indicators, one begins by transforming each date of 
occurrence of the annual floods to a Julian day (D;), that is, from O to 365; 
This implies not taking leap years into account. If a flood occurs on 
February 29, it is assigned the 28th. January dates remain the same, but 
February 31 is added, March 59, April 90, and so on until February. 
November which is added 304 and those of December are added 334, to 
obtain the Julian day. Next, the angle a, in radians corresponding to the 
date ¡ of each flood (D;) is obtained, with the following expression (Burn, 
1997; Cunderlik, Ouarda, 8: Bobée, 2004; Chen, Singh, Guo, Fang, € Liu, 
2013; Campos-Aranda, 2023b): 


a; =2m=X, with0 < a, < 21 (1) 


in which 
Tr = number pi with 3.141592654 as an approximate value 
X¡ = random variable of the dates of occurrence, in radians 


Next, the x and y coordinates of the flood occurrence dates 
described by the angles «a, are estimated based on the cosines and sines 


and their average values are obtained through the following equations: 


E= Ys cos (a;) (2) 
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y = 2YE, sen (ay) (3) 
where 


n = number of dates of occurrence of the annual floods analyzed 


Now, the mean direction (a) of the average date of the floods will 
be: 


x. 


a = arctan (2) (4) 


The application of the previous equation is done by first obtaining 
the arc tangent of y between x, both with a positive sign designated «, in 
radians; then if x and y are positive a=a, ifx < 0 and y > 0 x=u-a, if 
both are negative a =x+« and finally, if x > O and y < 0 a=2xr-a. The 
angles «; and « are converted to degrees (0% to 3609) by multiplying them 
by 57.295755. 


The value of « in Julian days is called mean flood day (DMC, by its 
acronym in Spanish), it is obtained by dividing by 27 and multiplying by 
365. The DMC index indicates the average occurrence date of the 
maximum annual flows in a given basin. Basins with similar DMC values 
can be expected to show similarities in other important hydrological 


characteristics. Logically, the DMC will be related to the size of the basin 


and its geographical location within the studied hydrological region (Burn, 
1997; Cunderlik et a/., 2004). 
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A measure of the variability of the n dates of occurrence of the 
floods, in relation to the DMC, can be estimated by calculating the 


resulting average, which is expressed as: 
r= IX +y (5) 


The seasonality index r is a dimensionless measure of the data 
dispersion, it takes values between zero and one. A unit value indicates 
that all the floods occur on the same date, while a value close to zero 


implies great variability of occurrences throughout the year. 


Ramirez-Orozco et al. (2009) establish the following five degrees of 
seasonality: (1) very strong, when 7 > 0.90, (2) strong, when r fluctuates 
between 0.70 and 0.90, (3) medium, when r varies from 0.50 to 0.70, 
(4) low, when r changes from 0.10 to 0.50 and (5) very low or weak, 
when r < 0.10. Chen et al. (2013) indicate that if r is close to unity, a 


single season or wet season can be expected to be dominant. 


The von Mises distribution 


This probabilistic model is commonly used to represent random variables 
that have two-dimensional direction and a single mode. For this reason, 
the von Mises distribution (dvM, by its acronym in Spanish) is considered 


the natural analogy of the Normal model for seasonal data. Its probability 
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density function is the following equation (Metcalfe, 1997; Carta, Bueno, 
8: Ramirez, 2008; Chen et a/., 2010): 


f() = EPlecosO a] with0<x<2r,0<u<2x,k>0 (6) 
2r Ip (xk) 


The dvM is symmetric with its mode at x = y, which is also ¡its mean 
direction (a) and the dispersion is given by the concentration parameter 
k (kappa). The denominator of Equation (6) makes the area under the 
curve unitary and for this reason it is called the normalization factor (FN, 
by its acronym in Spanish); includes the modified Bessel function of the 


first type of order zero [Zo(x)]. 


To estimate the non-exceedance probability of a value x, Equation 


(6) is numerically integrated, that is: 


FC) = 5d explx> cos(x =p] (7) 


2rrIp(x) 


The previous expression defines the PDF of the dvM, Jo(k) is 
estimated with the following ascending series, which comes from Olver 
(1972): 


2/4 12/47 2/4 2/4) 12/47 
E (5) 
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Numerical integration of the dvM 


To carry out such numerical integration of Equation (7), the Gauss- 
Legendre quadrature method was adopted, whose univariate operational 


equation is (Nieves 8 Domínguez, 1998; Campos-Aranda, 2003): 


SL F00 dx <p, wi: [Coorirore] A 


in which 
w¡= coefficients of the method whose abscissas are h; 


np = number of pairs in which the function f(x) is evaluated, with the 


argument indicated in f(:*) of Equation (9). 
In Davis and Polonsky (1972) the 12 used pairs of w; and h; with 15 


digits were obtained, because the Basic language accepts 16 digits as 


double precision variables. 


Fit of the dvM in the wet season 


When the annual flood occurrence dates cover, for the most part, a quite 
defined period in months, for example, from June to October, then the 
application of Equation (7) is carried out via numerical optimization, to 
find the values of y and k that reduce the sum of the differences between 


theoretical probabilities [F1(x)] and empirical probabilities [Fs*(x)] to the 


square, that ¡s: 
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Minimize FO = Y? ,[Fp(x) = F£(07? (10) 


The minimization of the previous objective function will be carried 
out by means of the Rosenbrock algorithm. The empirical probabilities are 
defined with the Gringorten formula (Chen et a/., 2010), which is: 


m-0.44 
n+0.12 


Fs(x) = 
(11) 


where 


m = order number of the data or occurrence date in radians (x = a¿), when 


they are located in progressive magnitude 
n = total number of data 


Logically, the occurrence dates outside the main period or wet 
season of the floods are previously eliminated, to improve the fit and, 
therefore, the definition of the dvM. More details of the previous process 


can be consulted in Campos-Aranda (2023b). 


Rosenbrock algorithm 


It is a numerical direct search procedure that attempts to define the 
minimum of a nonlinear function of multiple unbounded random variables. 


Rosenbrock's algorithrn assumes that the function is unimodal and begins 
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by defining a straight line, or search direction, from a given starting point. 
Next, it evaluates the objective function (FO) at various points on the line 
and determines the optimum; when this has happened, a new search 


direction is selected and the process is repeated recursively in stages. 


In this algorithm it is convenient to provide different starting points 
to search for the global minimum, based on the estimated local minima. 
A more detailed description of the process can be found in Rosenbrock 
(1960), Kuester and Mize (1973), and Campos-Aranda (2003). 


Advantages of Copula Functions 


As already indicated, the essential advantage of Copula Functions (FC, by 
its acronym in Spanish) consists in allowing to express a joint distribution 
of correlated random variables, as a function of their marginal 
distributions, previously adopted. So, a FC links or relates the univariate 
marginal distributions to form a multivariate distribution. Another basic 
advantage of FCs when forming multivariate distributions is the fact that 
they separate the effect of dependence between random variables from 


the effects of marginal distributions in joint modelling. 


Due to the above, the construction of the multivariate distribution 
is reduced to the study of the relationship between the correlated 
variables, if the univariate marginal distributions are known. The use of 
FCs offers complete freedom to adopt or select the univariate marginal 
distributions that best represent the data (Meylan et a/., 2012; Zhang 8 
Singh, 2019; Chowdhary 8 Singh, 2019). 
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Copula families to apply 


The Copula functions (FC, by its acronym in Spanish) that have been 
developed have been classified into four classes: Archimedean, extreme 
value, elliptic, and miscellaneous. They are also classified into one- 
parameter or multi-parameter copulas, depending on the amplitude with 
which the structure of the dependency between the variables X and Y is 
defined (Meylan et a/., 2012; Chowdhary 8 Singh, 2019). Salvadori et al. 
(2007) present a comprehensive and useful summary of FC that have 


been applied in the field of hydrology. 


Designating Fx(x) = u, Fy(y) = v and O the parameter that measures 
the dependence or association between u and v, we have the following 
two families of Archimedean copulas and extreme values (Salvadori et al., 
2007; Zhang € Singh, 2019; Chen 8 Guo, 2019; Chowdhary € Singh, 
2019). 


1. Frank. Its equation and variation space of O are: 


(e79U_4)(e7 811) 
e=0-4 


C(u,v) = Sn 14 (00, 00) Y (03 (12) 


For the negative dependence 0O<6<1 and for the positive one 0>1, with O 


= 1 for the independence between u and v. The relationship of O with r,, 


is the following: 
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Tn = 1+5[D,(0) — 1] (13) 
where D:(0) is the Debye function of order 1, expressed as: 


D,(0) =3 fas (14) 


es1 
The previous equation was estimated with numerical integration, 


ratifying its results with the values tabulated by Stegun (1972). 


2. Gumbel-Hougaard, which accepts only positive dependence. 


Its equation and variation space of O are: 
C(u,v) = exp(-[Elnw? + in ye, [1, 00) (15) 


With O = 1 there is independence between u and v. The relationship 


of 6 to Kendall's tau ratio is as follows: 


n= (16) 
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Association measures 
Concordance 


Since the FC characterizes the dependency between the random variables 
u and v, it is necessary to study the association measures, in order to 
have a method that allows estimating its parameter 6. In general terms, 
a random variable is concordant with another, when its large values are 
associated with the large values of the other and the small values of one 
with the small values of the other (Salvadori et al., 2007; Chowdhary € 
Singh, 2019). 


Some variables with direct linear correlation will be concordant, 
since as one increases the other does too. Variables with inverse linear 
correlation will be discordant, since large values of one will correspond to 
small values of the other. This implies that the pairs (x1-x2)(y1-y2)>0 are 
concordant (c) and discordant (d) when (x1-x2)(y1-y2)<0 (Salvadori et al., 
2007; Chowdhary €: Singh, 2019). 


A numerical measure of association is a statistic that indicates the 
degree of dependence or variable association. For comparison purposes, 
such measures range from zero to +1 or to -1, indicating perfect 
association positive at +1 or negative at -1. Kendall's tau ratio and 
Spearman's rho coefficient are two non-parametric measures that provide 
information about a special form of association or dependency, known as 
concordance (Salvadori et al., 2007; Chen 8 Guo, 2019). 
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Kendall's tau ratio 


It measures the probability of having matching pairs, which is why it is 
the quotient of c-d between c+d. The expression to estimate it with 
bivariate data is (Zhang €: Singh, 2006; Zhang € Singh, 2019): 


a 5 ua a sign|(x; — 3107 = y;)] (17) 


In the above equation: 

n = number of observations 

sign[*] = +1 if such pairs are concordant and -1 if they are discordant 
Genest and Favre (2007) present a test for the tau quotient, 

adopting Ho as the null hypothesis that X and Y are independent and then 

the statistic has an approximately Normal distribution with mean zero and 

variance 2(2n+5)/[9n(n-1)]. Then, Ho will be rejected with a confidence 


level a = 5 %if: 


9n(n—-1) 


ITnl > Zay/z = 1.96 (18) 


2(2n+5) 
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Dependency parameter estimation 


The simplest method to estimate the parameter 0 of the FCs is similar to 
the method of moments and is based on the inversion of the equation 
that relates O to the Kendall tau ratio or to Spearman's rho coefficient 
(Meylan et al., 2012; Chowdhary € Singh, 2019; Zhang € Singh, 2019; 
Chen 8 Guo, 2019). To obtain 6, in Equation (13) we proceed by trial and 


error; on the other hand, in Equation (16) its value is cleared. 


Estimation of joint empirical probabilities 


The bivariate empirical probabilities were estimated based on the 
Gringorten formula, applied by Yue (2000b), Zhang and Singh (2019), 
and Chen and Guo (2019). Such a formula is: 


j= ¡0.44 (19) 


— n+0.12 


in which: 
¡ = number of each piece of data, when they are ordered progressively 
n = total number of them, or the width of the processed record 


The previous expression was applied in the two-dimensional plane, 
with the data ordered progressively; the dates of occurrence (X;) in the 
rows and the maximum flows (Q;) in the columns. The plane formed is a 


square of n by n cells, with n cells on its main diagonal, when the order 
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number of the row is equal to that of the column. Then, each pair of 
annual data (X; and Q;) is located in the two-dimensional plane and the 
cell defined by the intersection of the row and column is identified with 


the number / that corresponds to the historical year drawn. 


When the n data pairs are drawn, year 1 is searched for and a 
rectangular or square area of minor X and Q values is defined, whose 
count of numbered cells within is NM, or minor X and Y combinations. 
Once the n values of NM, have been calculated, Gringorten's graphical 
position formula is applied to estimate the joint or bivariate empirical 


probability : 


Fay =PX<x0<q q == (20) 


n+0.12 


Selection of the Copula Function 


A simple approach to selecting the bivariate Copula function is based on 
the fit error statistics, by comparing the observed empirical probabilities 
(w?) with the calculated theoretical ones (wf) with the FC being tested. 
The indicators applied are the standard mean error (EME), the absolute 
mean error (EMA) and the maximum absolute error (EAM); their 


expressions are (Chowdhary € Singh, 2019): 


1 
n 


EME = |Ly% (w2-—w£) (21) 
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EMA = Yi l(w? —w£)| (22) 
EAM = max|(w — w¿)| (23) 


Dependency at the top end 


Generalities 


The most important criterion applied to select a bivariate FC is based on 
the magnitude of the dependence at the upper end of the joint 
distribution, which has an impact on the veracity of the extreme 
predictions. The upper-right-tail dependency (A,) is the conditional 
probability that Q is greater than a certain percentile(s) of Fo(q), given 
that X is greater than that percentile in Fx(x); as s approaches unity. The 
dependence on the lower-left-tail (A, ), compares Q to be less than X, when 


s approaches zero (Chowdhary € Singh, 2019; Salvadori et al., 2007). 


In relation to the exposed FCs, Frank Copulas have no tail 
dependencies. In contrast, the Gumbel-Hougaard copula has significant 


upper-tail dependency, equal to: 


ly = 2-28 (24) 


Dupuis (2007) tested six Copula families and found that their ability 


to estimate extreme events ranges from poor to good, in the following 
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order: Clayton, Frank, Normal, t-Student, Gumbel-Hougaard, and 
Asociated Clayton (Survival Clayton). Poulin, Huard, Favre and Pugin 
(2007) reached similar conclusions, when comparing the same six families 
of copulas and the so-called A12, which has significant right tail 


dependency. 


Estimation of the observed dependency 


In order to approach the estimation of the dependence in the upper tail 
(Au) shown by the available data, the so-called Empirical Copula must first 
be defined. Since the FC characterizes the dependence between the 
random variables X and Q, then the pair of ranks R; and S; coming from 
such variables is the statistic that retains the greatest amount of 
information and its scaling with the factor 1/(n+1) generates a series of 
points in the unit square [0,1]?, forming the domain of the Empirical 
Copula (Chowdhary 8 Singh, 2019), defined as follows: 


Cp (uv) = 2311 (4 <u LL <v) (25) 


n+1 7 O n+1 


In the above equation, 1(*) indicates a function of the random 
variables U and V, which are a continuously increasing transformation of 


X and Y, related to the empirical probability integrals Fr(X) and F,(Y), with 


the following equations: 
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U; = AS EL (xi) V¡= a Fr (q1) (26) 


n+1 n+1 


Poulin et al. (2007) use the estimator proposed by Frahm, Junker 
and Schmidt (2005), which is based on a random sample obtained from 


the empirical copula, its expression is: 


1 1 1 1 
12”. =2- Zexp (Ez in] Inn ¿In Gm (27) 


This estimator accepts that the FC can be approximated by one of 
the extreme values class, it has the advantage of not requiring a threshold 


value for its estimation. 


Ratification of the selected Copula function 


This is the most important stage of the FC application process, since it is 
verified that such model faithfully reproduces the observed joint 
probabilities (Equation (20)). Yue (2000a) describes a simple and 
practical way of representing the empirical and theoretical joint 
probabilities. It consists of taking the first one to the abscissa axis and 
the second one to the ordinate axis; in such a graph, each data pair 
defines a point that coincides with or departs from the 45% line. The 
inspection of the graph described and the value of the correlation 


coefficient —in these cases, greater than 0.98— ratify the validity of the 


joint probabilistic model. 
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Yue (2000b), and Yue and Rasmussen (2002) apply the 
Kolmogorov-Smirnov test with a significance level (a) of 5 %, to accept 
or reject the absolute maximum difference (dma, by its acronym in 
Spanish) between the joint probabilities. To evaluate the statistic (D,) of 
the test, the expression that Meylan et al. (2012), for a = 5 % this is: 


Da e 1.358 (28) 


n is the number of data. If the dma is less than D,, the adopted FC 


is ratified. 


Selection of marginal distributions 


The approach for selecting the marginal distributions was very simple and 
consisted of applying the three FDPs that have been established as 
reference or applicable under precept, which are the Log-Pearson type III 
(LP3), the General Extreme Values (GVE) and Generalized Logistics 
(LOG). In addition, three widely used models were applied: the 
Generalized Pareto (PAG), the Kappa and the Wakeby. The first four 


mentioned FDPs have three fit parameters and the last two, four and five. 


With the exception of the LP3 that was applied with the method of 
moments in the logarithmic (WRC, 1977) and real (Bobée, 1975) 
domains, the rest were adjusted with the method of L moments (Hosking 
es: Wallis, 1997). 
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The selection of the most convenient FDP was based on the value 
that each one generates with a probability of non-exceedance of 1 %; 
that is, a very low value that does not exceed the minimum observed 
values, to avoid negative marginal probabilities (vw). This selection 
criterion is required in records of maximum flow Q that present very low 


minimum values, compared to their maximum extremes. 


In addition, the standard fit errors (EEA) and absolute mean (EAM) 
were taken into account, as well as the magnitude of the predictions in 


the return periods greater than 500 years. 


Fitting errors 


The first criterion applied for the selection of the best FDP to some 
available data or series, were the so-called adjustment errors (Kite, 1977; 
Willmott 8 Matsuura, 2005; Chai € Draxler, 2014). This criterion and the 
one described to avoid negative probabilities will allow adopting an 
adequate distribution between the models: LP3, GVE, LOG, PAG, Kappa 
and Wakeby. 


By changing in equations (20) and (21), the probabilities observed 
by the ordered data of the analyzed series (x; or y;¡) and the probabilities 
calculated by the estimated values with the FDP that is tested or 
contrasted, the standard error of fit (EEA) and the mean absolute error 
(EAM) are obtained. The values that are estimated (%, or ),) are sought 
for the same probability of non-exceedance assigned to the data with the 


empirical Gringorten formula (Equation (19)). 
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Bivariate return periods 


The first bivariate return period of the event (X, Y) is defined under the 
OR condition, which indicates that the limits x or y, or both can be 
exceeded and then, the classical equation of the return period or inverse 
of the exceedance probability will be (Shiau, Wang, € Tsai, 2006; Genest 
8: Chebana, 2017): 


a e RNE (29) 


Txy = PX>AVY>y)  1-Fxy(0y)  1-C[FXGOFyO)] 


in which CLFXOO FyY(y)] is the selected FC. 


The second bivariate return period of the event (X, Y) is associated 
to the case in which both limits are exceeded (X>x,Y>y) or AND condition, 
its equation is (Shiau et al., 2006; Genest 8 Chebana , 2017): 


DS O ORAR (50) 


POSXAY>Y) Fx yy) 1 EXGOFyONACIEZO y ON] 


Aldama (2000) obtains the expression Fyy(x,y) of the bivariate 
probability of exceedance by means of a logical and simple probability 
reasoning applied in the Cartesian plane. Instead, Yue and Rasmussen 
(2002) resort to the Cartesian plane to define the bivariate event (X, Y) 


conceptually, which can occur in any of the four quadrants. 
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The relationship between bivariate and univariate return periods is 
as follows (Yue € Rasmussen, 2002; Shiau et al., 2006; Vogel € 
Castellarin, 2017): 


Txy < min[Tx, Ty] < max[Tx, Ty] < T'xy (31) 
being: 

1 1 
E Fx(x) — 1-Fx(x) (32) 


E e (33) 


Y HO) YO) 


The conditional bivariate return periods are based on the conditional 
distribution of X given that Y<y, which is expressed as follows (Chen 8 
Guo, 2019): 


PX < xlV < y) = CULO |FYO) < y) = EA (34) 


A similar equation is obtained for Y given that X<x. By similarity 


with the previous equation and Equation (30), the conditional distribution 


can be obtained for variables X and Y exceeding some limits (Chen € Guo, 
2019), this is: 
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_ P(X>x Y>y) ha 1-Fx CO -Fy(M+C[IFx CO Fy()] 
PA) (35) 


Data for processing 


In this study, the record of the Guamúchil hydrometric station, from 
Hydrological Region No. 10 (Sinaloa), Mexico, was processed. The 
Guamúchil gauging station has a basin area of 1 645 km? and a 32-year 
record that began in 1940 and ended in 1971, when the Eustaquio Buelna 
dam was built (Aldama et a/l., 2006). In this record, three years were 
eliminated: (1) 1949 whose flow of 375 m3/s occurred on January 22; (2) 
1960 with a rate of 1 373 m*/s occurred on January 13, and (3) 1968 with 
a rate of 200 m3/s that occurred on February 10. The record of the 


adopted wet season is shown in Table 1, in its columns 1 to 6. 


Table 1. Flows of annual floods, dates of occurrence, theoretical and 
empirical non-exceedance probabilities and count for joint empirical 
probability at the Guamúchil hydrometric station in Hydrological Region 


No. 10 (Sinaloa), Mexico. 


1 2 3 4 5 6 7 8 9 10 
i a; (radians) FIOO | Fe(x< 

No. | Y | mes [Día | DJ rod | FEO | 
(m3/s) obs ord ord ord 


1 205 AUG | 4 |216|3.718269 | 3.0641 | 0.041 |0.019| 3 


Z 65 SEP | 22 1265 |4.561765 | 3.2879 | 0.086 | 0.054 | 1 
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1 2 3]|415 6 7 8 9 [10 
i a; (radians) FIOO | Fe(x 
No. 2 Mes Día | DJ 0d FEO | y, 
(m3/s) obs ord ord ord 


3 | 445 |OCT| 2 |275|4.733907 3.5117| 0.165 | 0.088 | 12 
4 | 1550 | SEP | 26 |269|4.630622 | 3.5633 | 0.189 | 0.122 | 24 
[5 [| 392 |AUG| 30 |242|4.165838 | 3.5806 0.198 | 0.157 | 8 
[6 | 916 |OCT| 8 |281|4.837192|3.5806 0.198 | 0.191 | 22 | 
7 [| 241 |AUG| 10 |222|3.821554 | 3.6666 | 0.244 | 0.225 | 3 
(8 | 530 |AUG| 12 |224|3.855983 | 3.6838 0.254 | 0.260 | 9 | 
(9 [| 648 |JUL | 23 |204|3.511698|3.7183 0.275 |0.294 | 3 | 


10 272 AUG | 16 | 228 | 3.924839 | 3.7355 | 0.285 |0.328| 5 


11 422 SEP | 7 |250|4.303552 | 3.8216 | 0.342 [0.363 | 9 


12 377 AUG| 5 |217 3.735483 | 3.8560 | 0.366 | 0.397 | 4 


13 1179 SEP | 17 1260 4.475694 | 3.9076 0.402 | 0.431 | 21 


14 219 JUL | 10 [191 |3.287913 | 3.9248 | 0.415 [0.466 | 2 


15 3507 SEP | 23 |266 | 4.578979 | 3.9593 | 0.440 | 0.500 | 25 


16 165 JUN | 27 |178 3.064129 | 4.0281 | 0.491 [0.534 | 1 


17 526 AUG | 18 |230 | 3.959268 | 4.0970 | 0.543 | 0.569 | 10 


18 1014 SEP | 20 263 | 4.527337 | 4.1658 | 0.593 | 0.603 | 20 


19 1610 |AUG| 2 |214 3.683840 | 4.3036 0.689 | 0.637 | 8 


20 525 AUG | 1 |213|3.666626 | 4.3380 | 0.711 |0.672| 5 
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1 2 3 4 5 6 7 8 9 10 
i a; (radians) FIOO | Fe(x 
No. 2 Mes Día | DJ ro) [Febd |, 
(m3/s) obs ord ord ord 


21 | 985 [|OCT| 4 [277 4.768336 | 4.4413 | 0.772 0.706 | 22 
122 | 460 | JUL | 26 |207|3.563341 | 4.4757 | 0.791 | 0.740 | 3 
[23 | 390 |AGO| 26 |238|4.096982 | 4.5273 0.817 [0.775 | 7 
(24 | 449 | JUL | 27 |208|3.580555 | 4.5618 | 0.833 | 0.809 | 3 
[25 | 794 | JUL | 27 [208|3.580555 | 4.5790 | 0.841 | 0.843 | 6 
(26 | 720 |AUG| 22 |234|4.028124 | 4.6306 0.862 | 0.878 | 13 
127 | 312 |SEP| 9 |252|4.337981|4.7339 0.898 [0.912 | 6 


28 520 SEP | 15 |258 | 4.441266 | 4.7683 | 0.909 | 0.946 | 13 


29 1045 [|AUG| 15 | 227 | 3.907625 4.8372 | 0.927 [0.981 | 12 


DJ = Julian day. 

Obs = observed values. 

Ord = ordered values. 

FIRQOO = probability of theoretical non-exceedance (Equation (7)). 
FENO = probability of empirical non-exceedance (Equation (19)). 


NMi = no. of minor X; and Qi combinations (Equation (20)). 
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Wald-Wolfowitz test 


This non-parametric test has been exposed and applied by Bobée and 
Ashkar (1991), Rao and Hamed (2000), and Meylan et a/. (2012) to verify 
independence and stationarity in records of maximum annual flows (Q;). 
Therefore, it was proposed to apply this test to the records of a, in radians 


and of maximum annual flows, which must be random samples. 


Results and their discussion 
Finding marginal distributions 
Verification of randomness 


First, the randomness of the records to process was verified, based on the 
Wald-Wolfowitz Test, whose statistic (U) led to values of -0.720 and 
1.522, for the occurrence dates (a;,) in radians and for the flow maximums 


of Table 1. Since U is less than 1.96, both series or samples are random. 


Distribution of annual occurrence dates 


To apply the Rosenbrock algorithm, the data (columns 3 and 4) from Table 
1 of the Guamúchil station, were adopted as initial values of y and k, 4.25 
and 0.50, which define an initial FO of 1.040. After 15 stages and 83 
evaluations of the FO, the following was obtained: FO = 0.0456, uy = 
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4,04104, k = 3.79233, FN = 59.3770 and the results are concentrated in 
column 8 of Table 1 and Figure 1. This fit of the dvM covers the interval 
of occurrences from 3.0641 to 4.8372 radians, which correspond to the 
following dates: June 27 to October 8 (281 - 178 = 103 days). 


Probabilidad de no excedencia 
(adimensional) 


3.0 3.5 4.0 4.5 5.0 
Fechas de ocurrencia (radianes) 


Figure 1. Fit of the von Mises distribution to the annual flood 
occurrence dates recorded at the Guamúchil hydrometric station, in 
Hydrological Region No. 10 (Sinaloa), Mexico. X axis: Occurrence dates 


(radians). Y axis: Non—-exceedance probability (adimensional). 
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The estimation of dates corresponding to the return periods (Tr, 
years) that will contain the flood predictions, was explored by trial and 
error based on Equation (7), and the obtained values of y = 4.04104 and 
k = 3.79233 for recording Table 1. The month and day that define the 
most approximate values of non-exceedance probability [F(x)], related to 
the estimated Tr, are adopted. Table 2 shows the estimates obtained for 
the indicated Tr. 


Table 2. Estimated occurrence dates with von Mises distribution for the 
indicated return periods, in the annual floods recorded at the Guamúchil 


hydrometric station, Mexico. 


Tr (years) | Obtained Day in Estimate | Estimate Estimated 
assigned date radians d F(x) d 1-F(x) Tr (years) 
50 OCT 29 | 5.198691 | 0.97970 | 0.02030 49.3 
100 NOV8 | 5.370833 | 0.98969 | 0.01031 97.0 
500 NOV 28 | 5.715117 | 0.99806 | 0.00194 515.5 
1 000 DEC 4 | 5.818403 | 0.99904  0.00096 1041.7 


Distribution of maximum annual flows 


Table 3 shows the fit errors and predictions (m3/s) obtained with the three 
reference distributions and the three in general use, applied to the 
maximum flow records in Table 1. The minimum value of the maximum 


flows of 65 m3/s, was lower than the magnitudes with an exceedance 
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probability of 1 %, which are estimated with the FDP applied, except with 


the Kappa distribution, which was adopted accordingly. 


Table 3. Fit errors and predictions (m3/s) of the six distributions applied 
in the record of maximum annual flows of the floods in the Guamúchil 
hydrometric station, of the Hydrological Region No. 10 (Sinaloa), 


Mexico. 


EEA EAM Return periods in years 
FDP 


(m3/s) | (m?/s) 50 100 500 |1000 | 5000 | 10 000 


LP3 181.5 90.0 2649 | 3246 | 4855 5650 7747 8765 


GVE 152.6 FL 2749 | 3647 6852 8929 | 16361 | 21178 


LOG 140.4 72.2 2702 | 3649 7283 9798 | 19499 | 26219 


PAG 163,3 82.6 2738 | 3439 — 5493 6599 9837 11581 


KAP 2256 116.8 | 2430 | 3031 | 4946 6068 9664 11775 


WAK | 152.7 76.0 2789 | 3697 — 6839 8815 | 15651 | 19951 


The location (uz), scale (a2) and shape (k2 and h2) parameters of 
the adopted Kappa distribution are: 578.6213, 265.6867, -0.275 and - 
1.0, expressed in the following equation: 


1/h> 


F(y) = [12 [1 202012 (36) 
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The Kappa distribution adopts a value of minus one in its second 
shape parameter (h), because the L ratios (r¿ = 0.42744 and 71,= 0.32260) 
define in the L Ratio Diagram, a point above the Generalized Logistic 
distribution (Campos-Aranda, 2023a). 


Selection and ratification of the FC 


The bivariate data processing in Table 1 resulted in the following two 
association indicators: rxy = 0.3563 and r, = 0.2315. Equation (18), with 
n = 29 and the cited tau, gives a value of 0.9424; therefore, the Kendall 


ratio is not significant. 


Although the value of tau is not statistically different from zero, 
being positive indicates a direct correspondence or concordance, although 
low, between the flood occurrence dates and their maximum flow value, 
for their annual magnitudes. Such relationship or dependency will be 
modeled by the FC. 


Table 4 shows the statistical fit indicators that were obtained by 
applying the Frank and Gumbel-Hougaard (G-H) FCs. In equations (21) 
to (23), the empirical bivariate probabilities were estimated with Equation 
(20) and the theoretical ones with equations (12) and (15). 
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Table 4. Statistical indicators of the fit of the Copula Functions indicated 


in the annual floods of the Guamúchil hydrometric station, Mexico. 


No. | No. 
Copula 0 EME EAM MDP MDN Ay 
DP | DN 
Frank 2.1790 | 0.0309 | 0.0237 | 20 9 0.0682 | -0.0368 | 0.0000 
G-H 1.3013 | 0.0326 | 0.0246 | 18 11 |0.0785 | -0.0369 | 0.2965 


DP, DN = positive and negative differences. 


MDP, MDN = maximum positive and negative difference. 


Since the application of Equation (27) led to a AF“ value of 0.3416, 
there is no difficulty in selecting FC Gumbel-Hougaard in Table 4. 


On the other hand, since the value of 1£F% was slightly higher than 
the 2, of the FC, one should look for a FC with greater dependence on its 
right end, which was not necessary in this case, because predictions were 
not made for high return periods, but rather estimates of exceedance 
probabilities were formulated, according to occurrence dates, based on 
conditional return periods. Chen and Guo (2019) exclusively apply the 


Gumbel-Hougaard FC, in this type of frequency analysis. 


Table 5 shows the observed empirical bivariate non-exceedance 
probabilities (w?) calculated with Equation (20) and estimated theoretical 


(w?) with the Gumbel-Hougaard FC. The maximum positive and negative 


differences are also indicated shaded. 
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Table 5. Joint non-exceedance probabilities and their differences, 


calculated with the Gumbel-Hougaard FC, for the annual floods of the 


Guamúchil station, Mexico. 


No. wW; w; Differences 
1 0.0879 | 0.0785 0.0094 
2 0.0192 | 0.0561 -0.0369 
E 0.3970 | 0.3529 0.0441 
4 0.8091 | 0.8232 -0.0141 
5 0.2596 | 0.2318 0.0279 
6 0.7404 | 0.7201 0.0203 
7 0.0879 | 0.0881 -0.0002 
8 0.2940 | 0.2155 0.0785 
9 0.0879 | 0.1215 -0.0335 
10 | 0.1566 | 0.1165 0.0400 
11 0.2940 | 0.2812 0.0128 
12 0.1223 |.0,1229 -0.0007 
13 | 0.7060 | 0.7126 -0.0066 
14 | 0.0536 | 0.0250 0.0285 
15 | 0.8434 | 0.8391 0.0043 
16 | 0.0192 | 0.0103 0.0090 
17 0.3283 | 0.2514 0.0769 
18 | 0.6717 | 0.6919 -0.0202 
19 | 0.2596 | 0.2485 0.0111 
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No. wi wj Differences 


20 | 0.1566 | 0.1486 0.0080 


21 | 0.7404 | 0.7410 -0.0007 


22 | 0.0879 | 0.1043 -0.0163 


Za | 0.2253 | 0.2156 0.0097 


24 | 0.0879 | 0.1061 -0.0182 


25 | 0.1909 | 0.1628 0.0282 


26 | 0.4313 | 0.3619 0.0694 


27 | 0.1909 | 0.1998 -0.0089 


28 | 0.4313 | 0.3865 0.0448 


29 1 0.3970 1:0,3623 0.0346 


On the other hand, Equation (28) defines D, = 0.2522 and since the 
maximum absolute difference in Table 5 is 0.0785, the Kolmogorov- 
Smirnov test ratifies the adopted Gumbel-Hougaard (G-H) FC. The 
correlation coefficient (rxy) between the empirical and theoretical 
probabilities, estimated with the FC of G-H, was 0.9931; therefore, the fit 
is excellent. The graphic contrast between both probabilities, to ratify the 
adoption of FC of G-H, is shown in Figure 2 for the data in Table 5. 
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Figure 2. Graphic contrast of joint probabilities of the occurrence dates 
and the maximum flows of the floods recorded at the Guamúchil 
hydrometric station, Mexico; with the Gumbel-Hougaard FC. X axis: 
Empirical non-exceedance probability. Y axis: Theoretical non- 


exceedance probability. 
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Graphs of joint return periods 


The bivariate return periods of the AND type are calculated with Equation 
(30). For their estimates with values of Tyy of 50, 100 and 500 years, a 
date (month and day) is defined and its non-exceedance probability F(x) 
with Equation (7) and the values of y = 4.04104 and k = 3.79233. Now 
by trial and error, a maximum flow (y = Qmax) is assigned and its F(y) is 
estimated with Equation (36). Both results are taken to Equation (15), 
with O = 1.3013, to estimate the joint non-exceedance probability and 
thus obtaining with Equation (30), the 7T;, that must coincide with the 


searched value. 


Dates of occurrence and maximum flows are selected arbitrarily, in 
order to define the curves of Txy. Table 6 shows the values adopted to 


define the three graphs in Figure 3. 


Table 6. Pairs of dates of occurrence and maximum annual flow used to 
define the graphs of the AND type joint return period with the Gumbel- 


Hougaard FC, in the floods of the Guamúchil station, Mexico. 


Txy = 50 years Txy = 100 years Txy = 500 years 
m m m 

Date ND* Q Date ND* Q Date ND* Q 
M3/s M3/s M3/s 


Oct 29 121 117 Nov 8 131 255 Nov 28 151 0 


Oct 27 119 569 Nov 6 129 717 Nov 27 150 607 


Oct 25 117 814 Nov 5 128 | 880 Nov 26 149 | 1055 
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Txy = 50 years Txy = 100 years Txy = 500 years 
m m m 
Date ND* Q Date ND* Q Date ND* Q 
M3/s M3/s M3/s 


Oct 22 114 1097 Nov 3 126 | 1159 | Nov 25 148 1452 


Oct 20 112 1248 Oct 31 123 | 1496 | Nov 24 147 1798 


Oct 15 107 1530 Oct 28 | 120 | 1749 | Nov 23 146 | 2095 


Oct 10 102 1719 Oct 25 | 117 | 1941 | Nov 22 145 | 2347 


Oct 5 97 1854 Oct 20 | 112 | 2171 | Nov 21 144 | 2560 


Oct 1 93 1937 Octi5 | 107 | 2331 | Nov 20 143 | 2743 


Sep 25 87 2034 Octi0 | 102 | 2450 | Nov 18 141 | 3036 


Sep 20 82 2098 Oct 5 97 2543 | Nov 15 138 | 3357 


Sep 15 77 2152 Sep 30 92 2617 | Nov 10 133 | 3707 


Sep 15 72 2197 Sep 25 87 2679 Nov 5 128 | 3937 


Sep 5 67 2235 Sep 20 82 2731 | Oct 31 123 | 4103 


Aug 31 62 2268 Sep 15 77 2776 | Oct 25 117 | 4252 


Aug 25 56 2302 Sep 10 72 2815 | Oct 20 112 | 4349 


Aug 20 al 2326 Sep 5 67 2849 | Oct 15 107 | 4429 


Aug 15 46 2346 Aug 31 62 2879 | Oct 10 102 | 4496 


Aug 10 41 2363 Aug 25 56 2910 Oct 5 97 4554 


Aug 5 36 2317 Aug 20 51 2932 | Sep 30 92 4604 


Jul 31 31 2389 Aug 15 46 2951 | Sep 25 87 4648 


Jul 25 25 2400 Aug 10 41 2967 | Sep 20 82 4687 


Jul 13 13 2415 Aug 5 36 2980 | Sep 15 77 4722 


Tecnología y ciencias del agua, ISSN 2007-2422, 


Open Access bajo la licencia CC BY-NC-SA 4.0 (https: //creativecom- 15(4), 80-136. DOI: 10.24850/j-tyca-2024-04-03 
mons.org/licenses/by-nc-sa/4.0/) 4 : : : yy 


a 1) Check for updates 
OPEN ACCESS 


Tecnología y 
CienciaszAgua 
Txy = 50 years Txy = 100 years Txy = 500 years 
m m m 
Date ND* Q Date ND* Q Date ND* Q 
M3/s M3/s M3/s 
Jul 1 1 2430 Jul 31 31 2992 Sep 10 72 4754 


Jul 25 25 3002 Sep 5 67 4783 


Jul 20 20 3009 | Aug 31 62 4808 


Jul 15 15 3015 | Aug 20 5d 4855 


Jul 10 10 3019 | Aug 10 41 4886 


Jul 1 nl 3031 Aug 5 36 4899 


Jul 31 E! 4909 


Jul 25 25 4919 


Jul 20 20 4926 


Jul 1 1 4946 


*ND = day number from the 1st. of July. 
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Fechas de ocurrencia en número de día 
Figure 3. Graphs of the three design joint return periods Ty, obtained 
with the Gumbel-Hougaard FC, in the floods of the Guamúchil 
hydrometric station, Mexico. X axis: dates of occurrence in day number. 


Y axis: maximum annual flow (m/s). 


Probabilities of occurrence of design events 


Table 7, Table 8 and Table 9 show the conditional exceedance probabilities 
P(X>x|Qm>q0), calculated with Equation (35), when qo has a return 
period (Tr) of 50, 100 and 500 years, that is, when the marginal 
probability v = Fr(y) is equal to 0.98. 0.99 and 0.998 and the flow Qm 
exceeds 2 430, 3 031 and 4 946 m/s, according to Table 3. 
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dates provided that the maximum flow has a Tr = 50 years in the floods 


of the Guamúchil hydrometric station, Mexico. 


1 2 3 4 5 6 7 8 9 10 11 12 13 14 
Proposal PE (%) 99 95 90 80 70 60 50 40 30 20 10 5 1 
Date obtained JUN 3 | JUN 30 [JUL 13|JUL 27| AUG 6 | AUG 15 | AUG 23 | AUG 30 | SEP 8 [SEP 18| OCT 2 | OCT 15 | NOV 8 
a (radians) = X; |2.6510| 3.1158 |3.3396|3.5806|3.7527| 3.9076 | 4.0453 | 4.1658 |4.3208 |4.4929|4.7339| 4.9577 |5.3708 
Real PE 0.9902| 0.9509 |0.8994|0.8025|0.7038| 0.5977 | 0.4959 | 0.4068 |0.3000|0.2002|0.1016| 0.0488 |0.0103 
PNE of the FC GH |0.0097| 0.0489 |0.1002|0.1967|0.2949| 0.4003 | 0.5014 | 0.5897 |0.6954|0.7939|0.8902| 0.9404 |0.9739 
Conditional PE 0.9985| 0.9916 |0.9813|0.9592|0.9333| 0.9012 | 0.8652 | 0.8280 |0.7726|0.7016|0.5853| 0.4591 |0.2106 

PE = probability of exceedance, dimensionless. 
PNE = probability of non-exceedance, dimensionless. 
GH = Gumbel-Hougaard FC. 
Table 8. Calculations of the conditional probability of the occurrence 
dates provided that the maximum flow has a Tr = 100 years in the 
floods of the Guamúchil hydrometric station, Mexico. 

1 2 3 4 5 6 7 8 9 10 11 12 13 14 
Proposal PE (%) 99 95 90 80 70 60 50 40 30 20 10 5 1 

Date obtained JUN 3 | JUN 30 [JUL 13 |JUL 27| AUG 6 | AUG 15 | AUG 23 | AUG 30 | SEP 8 [SEP 18| OCT 2 | OCT 15 | NOV 8 
a (radians) = X; |2.6510| 3.1158 |3.3396|3.5806|3.7527| 3.9076 | 4.0453 | 4.1658 |4.3208 |4.4929 |4.7339| 4.9577 |5.3708 
Real PE 0.9902| 0.9509 |0.8994|0.8025|0.7038| 0.5977 | 0.4959 | 0.4068 |0.3000|0.2002|0.1016| 0.0488 |0.0103 
PNE of the FC GH |0.0098| 0.0490 |0.1005|0.1972|0.2957| 0.4015 | 0.5030 | 0.5918 |0.6982|0.7974|0.8951| 0.9467 |0.9828 
Conditional PE [|0.9988| 0.9932 |0.9849|0.9671|0.9461| 0.9202 | 0.8911 0.8609 |0.8160/|0.7582|0.6622| 0.5543 |0.3073 


PE = probability of exceedance, dimensionless. 


PNE = probability of non-exceedance, dimensionless. 
GH = Gumbel-Hougaard FC. 
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Table 9. Calculations of the conditional probability of the occurrence 
dates provided that the maximum flow has a Tr = 500 years in the 


floods of the Guamuchil hydrometric station, Mexico. 


1 2 3 4 5 6 7 8 9 10 11 12 13 


Proposal PE (%) 99 95 90 80 70 60 50 40 30 20 10 5 


Date obtained JUN 3 | JUN 30 [JUL 13 | JUL 27 | AUG 6 | AUG 15 | AUG 23 | AUG 30 | SEP 8 [SEP 18| OCT 2 | OCT 15 


a; (radians) = X; |[2.6510| 3.1158 |3.3396|3.5806|3.7527| 3.9076 | 4.0453 | 4.1658 |4.3208 | 4.4929|4.7339| 4.9577 


Real PE 0.9902| 0.9509 |0.8994|0.8025/0.7038| 0.5977 | 0.4959 | 0.4068 |0.3000|0.2002|0.1016| 0.0488 


PNE of the FC GH |0.0098| 0.0491 |0.1006|0.1975|0.2961| 0.4022 | 0.5039 | 0.5930 |0.6997|0.7995|0.8980| 0.9506 


Conditional PE |0.9993| 0.9958 [|0.9907|0.9798 0.9670| 0.9510 | 0.9332 | 0.9147 |0.8870|0.8514|0.7919| 0.7231 


PE = probability of exceedance, dimensionless. 
PNE = probability of non-exceedance, dimensionless. 
GH = Gumbel-Hougaard FC. 


Then, for the case of go = 2 430 m/s (Table 7), the following event 
probabilities can occur: (1) the occurrence probability of such flows after 
July 13 is 98.13 %; (2) the probability that such flows occur between July 
13 and November 8 will be: 98.13 - 21.06 = 77.07 %, and (3) during the 
period from September 8 to 18 will be: 77.26 - 70.16 = 7.10 %. 


Then, for the case of go = 3 031 m?/s (Table 8), the following event 
probabilities can occur: (1) the occurrence probability of such flows after 
July 13 is 98.49 %; (2) the probability that such flows occur between July 
13 and November 8 will be: 98.49 - 30.73 = 67.76 %, and (3) during the 
period from September 8 to 18 will be: 81.60 - 75.82 = 5.78 %. 


Finally, for the case of go = 4946 mY/s (Table 9), there are three 
following event probabilities (1) the probability of occurrence of such flows 


after July 13 is 99.07 %; (2) the probability that such flows occur between 
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July 13 and November 8 will be: 99.07 - 54,23 = 44.84 %, and (3) during 
the period from September 8 to 18 will be: 88.70 - 85.14 = 3.56 %. 


The comparison between the results of the three limits established 
for the maximum flow (qo), indicates that the probability of occurrence 
after July 13 increases slightly as the qo grows, which is logical, given the 
positive correlation of rxy = 0.3563 shown between the occurrence dates 
and the maximum flows. The above is also influenced by the denominator 


of Equation (35); which, as qo increases, is reduced from 0.02 to 0.002. 


On the contrary, the two occurrence probabilities calculated for the 
two periods within the wet season decrease as go increases, which is due 
to the lower probability of exceedance of each maximum flow, having 2, 
1 and 0.2 %. 


Contrast with results of a previous work 


Chen and Guo (2019) present on their pages 42 to 44 and Table 3.2, as 
a Case study of a bivariate flood frequency analysis of occurrence dates 
and maximum flow, its application in the multi-purpose Geheyan 
reservoir, with a basin of 17 000 km?. which receives an average annual 
rainfall of 1 500 millimeters and whose flood season covers five months 
from May 1 to September 30 (153 days). Keeping the due proportions, 
the results of Table 3.2 coincide with those of Table 8, since both are for 


a maximum flow of a return period of 100 years. 
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Application in basins with two wet seasons 


In general, in large basins, only the storms which result in cold fronts and 
hurricanes, are covered entirely with low intensity and long duration rains, 
which generate floods of great magnitude, due to their maximum flow, 
volume and duration. On the other hand, in such basins, convective 


storms are local, of great intensity and give rise to ordinary floods. 


In Mexico, the large basins of the mountainous zone of Hydrological 
Region No. 10 (Sinaloa) present two dates of occurrence of their floods, 
in autumn and winter (November to March) and in summer (June to 
September). In these flood regimes there are two modes and their 
probabilistic characterization of their occurrence dates is carried out with 
a mixture of von Mises distributions (Carta et a/., 2008; Campos-Aranda, 
2023b). 


Chen et al. (2010), and Chen and Guo (2019) have presented floods 
seasonal analyzes to define times of occurrence, in basins with two or 


more wet seasons. 


Application in large basins with a wet season 


In large basins with flood regimes of a single wet season, the bivariate 
study of dates of occurrence and maximum flow, in their various 
tributaries or main rivers can help to understand the evolution or 


development of their floods and then make more reliable forecasts of 


arrival dates downstream. 
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Conclusions 


The bivariate frequency analysis of dates of occurrence (X) and maximum 
flow (Qm) of annual floods is feasible due to the Copula functions; which 
allow building their joint probability distribution, based on some univariate 
marginals. The von Mises distribution is the one that characterizes the 


occurrence dates, and the Om in particular is an ideal FDP. 


In this study, a copula function (FC) of a single fitting parameter (0) 
was used, which is estimated based on the Kendall tau ratio, which is 
calculated with the joint registration of X and Qm. Such an approach first 
estimates AF“, or observed in the right tail dependency of the available 
record set. Afterwards, an FC is that reproduces such a value of A$"“ is 
searched; in this study the Gumbel-Hougaard was chosen. An FC that 
does not have a significant right tail dependency, such as Frank's, is also 
applied to allow comparing and assessing the quality of the previously 
adopted FC fit. 


The described numerical application, in the 32 annual data of 
occurrence dates and maximum flow (m3/s) of the annual floods recorded 
in the Guamúchil hydrometric station of the Hydrological Region No. 10 
(Sinaloa), Mexico; showed in Figure 2 a reliable reproduction of the 
empirical and theoretical bivariate probabilities, through the Gumbel- 


Hougaard FC, with a linear correlation coefficient of 0.9931. 


On the other hand, in Figure 3, with regard to the joint return 
periods of AND type design, infinite pairs of critical X and Qm can be 


defined, since they are in the curved region of each graph. 
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Finally, Tables 7, 8 and 9 show the calculation of conditional 
probabilities, which allow estimating, for an arbitrarily adopted period, the 
exceedance probability of design events with return periods of 50, 100 


and 500 years. 
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